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CROSS REFERENCE TO RELATED APPLICATIONS 



[001] This application claims priority from United States Provisional App. Ser. No. 
60/416,068 filed on October 4, 2002. 

5 

FIELD OF THE INVENTION 



[002] This invention relates to the field of seismic data processing and, more 
particularly, to creating velocity models for use in seismic data migration for determining 
10 subsurface earth structure represented by a 3-D volume of data for identifying structural 
and stratigraphic features in three dimensions. 



BACKGROUND OF THE INVENTION 

1 5 [003] Searching for subsurface mineral and hydrocarbon deposits comprises data 

acquisition, analysis, and interpretation procedures. Data acquisition involves energy 
sources generating signals propagating into the earth and reflecting from subsurface 
geologic structures. The signals received are recorded by receivers on or near the surface 
of the earth. The received signals are stored as time series (seismic traces) that consist of 

20 amplitudes of acoustic energy which vary as a function of time, receiver position, and 

source position and, most importantly, vary as a function of the physical properties of the 
structures from which the signals reflect. The data are generally processed to create 
volumes of acoustic images from which data analysts (interpreters) create maps and 
images of the subsurface of the earth. 
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[004] Data processing involves procedures that vary depending on the nature of the data 
acquired and the geological structure being investigated. A typical seismic data 
processing effort produces images of geologic structure. The final product of data 
5 processing sequence depends on the accuracy of these analysis procedures. 

[005] Processed seismic data are interpreted to make maps of subsurface geologic 
structure to aid decisions for subsurface mineral exploration. The interpreter's task is to 
assess the likelihood that subsurface hydrocarbon deposits are present. The assessment 
10 will lead to an understanding of the regional subsurface geology, important main 

structural features, faults, synclines and anticlines. Maps and models of the subsurface, 
both in 2D and 3D representations are developed from the seismic data interpretations. 
As is well known in the art, the quality and accuracy of the seismic data processing has a 
significant impact on the accuracy and usefulness of the interpreted data. 

15 

[006] High quality data processing greatly simplifies data interpretation, since resources 
can be focused on the geologic structure since subsurface imaging can be made less 
ambiguous. Unfortunately, three dimensional geophysical data processing and/or 
modeling frequently require large computation expenses, and practitioners are forced to 
20 simplify the data processing effort as much as possible to reduce analysis time and cost. 

[007] The sheer volume of data impacts data processing considerations. Seismic survey 
data sets can involve hundreds of thousands of source locations, with each source 
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location associated with many hundreds more receiver locations. Each input/output data 
transfer demand burdens resources independent of the computation burden. 

[008] There have been several different approaches to manage these computational 
5 resource burdens. These approaches relate to the manner in which the data acquisition 
exercise is designed and carried out, as well as to assumptions made during data 
processing. The use of available a priori geologic and geophysical information can 
facilitate the minimization of the seismic data acquisition effort. Such a minimization of 
resources reduces the amount of data that is acquired by reducing the acquisition effort. 

10 

[009] Minimization of the computational effort is often implemented during data 
processing. Compromises made during data acquisition and/or processing may lead to 
ambiguous and/or inaccurate images. Because little is generally known of the geologic 
structure being investigated, the interpreter will not know the extent to images are 
15 erroneous. 

[0010] It is not uncommon for significant computer resources to be involved when large 
or complex data volumes are processed, often involving weeks or months of actual 
computer processing time. The recent availability of massively parallel processor 
20 computers offers a significant opportunity to reduce overall processing times. Massively 
parallel processors (MPPs) can have multiple central processing units (CPUs) which can 
perform simultaneous computations. By efficient use of these CPUs, projects that took 
weeks or months of resource time previously can be reduced to a few days or a few 
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hours. These advantages can be enhanced further when efficient algorithms are included 
in the MPP software. 

[0011] Computational algorithms have previously been written for prior seismic analysis 
5 routines using single or just a few processors, usually using sequential computing. 

Sequential computing performs single procedures at any given time. Options for 
obtaining enhanced performance are limited when few processors are available. 

[0012] MPP computing machines offer an obvious computation advantages. The total 
10 time required to process a dataset can be reduced by dividing the work to be done among 
the various CPUs or CPU clusters in manner such that each CPU performs useful work 
while other CPUs also work in parallel. 

[0013] From Robert E. Sheriffs "Encyclopedic Dictionary of Exploration Geophysics" 
15 SEG Press, 2002: Tomography is a method for finding the velocity and reflectivity 

distribution from a multitude of observations using combinations of source and receiver 
locations, or of determining the resistivity distribution from conductivity measurements 
using a transmitter in one well and a receiver in another well. Tomography is derived 
from the Greek for "section drawing." Generally space is divided into cells and the data 
20 are expressed as line integrals along raypaths through the cells. Transmission tomography 
involves borehole-to-borehole, surface-to-borehole, or surfaceto-surface observations. 
Reflection tomography involves surface-to-surface observations as in conventional 
reflection or refraction work. In seismic tomography, slowness or velocity, and 
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sometimes an attenuation factor, is assigned to each cell and traveltimes and amplitudes 
are calculated by tracing rays through the model. The results are compared with observed 
times and amplitudes; the model is then perturbed and the process repeated iteratively to 
minimize errors. Raypaths have to be recalculated after each change of assumed velocity. 

5 Diffraction tomography involves calculations assuming least-time travelpaths according 
to Fermat's principle rather than Snell's law bending at cell boundaries. Layer-based 
tomography divides the earth into layers, allowing for lateral variation of velocity within 
the layers, instead of subdivision into cells. Tomographic methods include the algebraic 
reconstruction technique, the simultaneous reconstruction technique, and Gauss Seidel 

10 methods. 

[0014] Tomography is used to compute corrections to velocities from observed 
traveltime errors in seismic datasets. The position of events observed on post-migration 
common image point gathers indicates errors that can be caused by several reasons: error 

15 in velocity, error in the depth location of reflector, or error in placing a fault. In some 
cases the velocity will be over corrected by assuming all the observed traveltime errors 
originated by velocities. Over corrected velocities will force more structural errors in 
deeper reflectors, and errors repeat and can increase as depth increases so that a correct 
velocity model cannot be obtained from further iteration. If the errors caused by 

20 misplaced horizons and fault locations can be corrected, the result is a better depth image 
and a better correction to the velocity model. 



COR-1060 



6 



[0015] For large datasets parallel computation methods offer significant processing time 
savings, and this is true for velocity analyses as well. An option for overcoming increased 
computational expenses is to employ efficient parallel algorithms and techniques. 

5 [0016] It would therefore be desirable to have a system and method that is able to 

perform tomographic velocity analyses efficiently on a parallel computer in a cost- 
effective manner. The present invention satisfies this need. 

SUMMARY OF THE INVENTION 

10 

[0017] The present invention provides a method and system for distributed residual 
tomographic velocity analysis using dense residual depth difference maps. Prestack 
seismic imaging is performed using an initial velocity field and interpreted horizons. A 
residual depth difference is estimated referenced to fixed offset and all horizons. 

15 Residual depth difference maps (p-maps) are computed for each offset and each horizon. 

The residual depth difference maps are back projected to determine slowness 
perturbation. The initial velocity model may be converted to slowness and the estimated 
slowness is composited therewith to produce a new slowness volume. The new slowness 
volume is converted to a new velocity volume for performing prestack seismic imaging. 

20 This process is repeated until the slowness perturbation is negligible or reaches a 
predetermined threshold. 
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BRIEF DESCRIPTION OF THE DRAWINGS 



[0018] The present invention and its advantages will be better understood by referring to 
the following detailed description and the attached drawings in which: 

Figure 1: Illustrates a flow chart of incremental velocity updating. 

Figure 2: Illustrates model partitioning and model sharing. 

Figure 3: Illustrates the general structure of distributed tomography. 

Figure 4 A: illustrates a flow chart of an embodiment of the invention. 

Figure 4B: illustrates a flow chart of another embodiment of the invention. 

Figure 5: illustrates a computer system for carrying out an embodiment of the invention. 

[0019] While the invention will be described in connection with its preferred 
embodiments, it will be understood that the invention is not limited thereto. On the 
contrary, it is intended to cover all alternatives, modifications, and equivalents which 
may be included within the spirit and scope of the invention, as defined by the appended 
claims. 

DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS 

[0020] A good macro velocity model is the key to producing good depth migration 
results. When the velocity is correct, Common Image Gather (CIG) reflections from 
different offsets should be aligned at the same depth. While this disclosure explains the 
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concepts using CIGs, the method and system of the present invention is directly 
applicable to Common Angle Image Gathers (CAIG) as well When the velocity is 
incorrect, it can be updated using the depth error in the CIG. 

[0021] 3D velocity model updating typically requires several gigabytes of memory. This 
is usually beyond the memory limit of a single modern computer. On PC clusters, large 
problems can be partitioned into small ones. Each partition can be processed 
independently. The final velocity is the combination of all the local models. 

[0022] Distributed tomography on a PC cluster is very cost effective. It makes it possible 
to solve extremely large problems quickly. In addition, very large offsets can be used in 
the inversion process. 

[0023] A macro velocity model is needed for depth migration. The kinematic information 
provided by this smooth velocity is the traveltime in seismic data and depth information 
in depth images. Depth migration and velocity estimation are coupled problems. Given a 
good velocity model, modern depth imaging algorithms can give satisfying results. 
Velocity estimation or inversion has been one of the most challenging geophysical 
problems. There are many ways to estimate a smooth background velocity. Prestack 
migration is one of them. 

[0024] One approach to velocity updating is described in Bednar et al., 1999, where the 
depth image is sorted into CIG gathers in which at each point there is a three-parameter 
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function: depth, offset and CDP. The principle is that reflections from one reflector at 
different offsets in each CDP should be mapped into one horizontal event in the depth 
image if there is no error in the input velocity model. For real problems, velocity models 
are never exact. The error in velocity usually appears as curved events in the CIGs. To 
determine the velocity, it is necessary to invert the depth error of each of these points in 
the CIG gathers. 

[0025] Ray-tracing is use to compute the difference between the input traveltime and 
traveltime computed form the model. This residual is used to update the velocity model. 
Information of each ray segment needs to be recorded. Obviously, for typical 3D real 
problems, a large amount of memory must be used. It is almost impossible to solve the 
typical problems on a single computer node. 

[0026] Herein disclosed is an implementation of 3D residual tomography on a distributed 
PC cluster. The full problem is divided into small pieces or partitions. Each small 
partition is solved on different node. The main issues in this method are model partition 
and consistency be tween overlapping pieces of the desired model. 

[0027] Residual Tomography CIG gather-based residual tomography updates the velocity 
model incrementally. The process is illustrated in the flow chart of Figure 1. An initial 
velocity 101 is used to migrate 103 seismic data to get a common image gather 105. A 
determination is made whether the CIG is flattened 107. If the CIG is not sufficiently 
flattened tomography 109 is performed to derive a new velocity 111 for input to the 



COR-1060 



10 



migration 103. The cycle is repeated until the CIG is sufficiently flattened 107 and the 
process finishes 113. 



10 



15 



[0028] For each iteration, the following optimization problem is solved 
C(m) = £^|4m-af (1) 



i=l 



[0029] Here, A ( is the linear operator that transforms the model m to the data space. 
Parameter bi is the known data. Parameter is a damping factor to ensure the 
smoothness and stability of the solution of the linear system below. Solving Equation 1 
is equivalent to solving the following linear system 
A T Am = A T b (2) 
where 



A = 



and 



MA 
M 2 A 2 

ju>A 2 



Iterative methods such as conjugate gradient (Paige and Saunders, 1982) can be used to 
solve this large linear system. 



[0030] Distributed Tomography: Residual tomography is highly automatic. Given 
model and data, the process can proceed until a solution is reached. The major problem 
20 with residual tomography is that input and output models are usually so large that they 
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can not fit into the memory of one computer node. In addition to memory resources 
required for the model, ray tracing, used for traveltime computation, also requires a large 
amount of memory. In 3D the following parameters are stored: the coordinates of the two 
endpoints of a ray segment and the direction vector of a ray segment. The ray tracing is 
5 done for each depth point of each offset of each CDP, for each inline and for each 
crossline. 

[0031] This is a five dimensional data space. To reduce the computation, one usually 
limits the number of xlines and inlines or limits the range of offsets. But too small ranges 
10 of offset or too few lines will not give reasonable or good velocity estimates. The best 
way is to divide the whole model into small partitions, each of which is solved in one 
computer or node. 

[0032] The reason we can do this is that, based on the locality principle, one point in the 
15 model space only influences on the points which fall in a circle on the surface. The radius 
of the circle is half of the largest offset for that CDP. Thus, we only need to consider the 
data within that circle in the velocity inversion of that point. This is illustrated with the 
portioned survey 201 of Figure 2. For each simple partition 203, additional areas are 
padded in all the four sides in (x,y) domain as illustrated by the offset buffer zone 205 in 
20 Figure 2 that represents half the offset of a seismic gather. These padded areas 205 are to 
be shared by the adjacent nodes. 
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[0033] Model Space Partition and Model Sharing: A key problem with distributed 
tomography is how to partition the model. The partitioning method depends on the 
characteristics of the model and the input depth residuals. Once partitioned, on each node, 
the space used by the model is relatively small compared with the space used by the :ay 
5 segments. Thus, the ideal partition can be based on the density of the surface depth 

residual estimates. The higher the density, the smaller the model partition. This usually 
results in load balanced computation. In addition to this complicated partition method, 
other simpler methods can be used. For example, if the model is much longer in inline 
direction than in crossline direction, we can partition the model in inline direction. 

10 

[0034] The partition used here is an intelligent one. The user specifies the number of 
partitions or the number of nodes, this number is factorized into two factors which are as 
close as possible. For example, 12 is factorized into 3 and 4 instead of 2 and 6. Then the 
larger number of these two factors is assigned to the longer direction. 

15 

[0035] Simply dividing the entire model evenly and assign each piece to each node won't 
work well because the part of the model on one node may need data which resides on an 
adjacent node. In order to make sure that all the data needed for the inversion of every 
point in each partition are used, each piece needs to be padded before inversion. This is 
20 shown in Figure 2. Partitions after padding are illustrated 207. 

[0036] Velocity Smoothing: Smooth velocity models are critical for obtaining good 
depth images, especially for prestack depth migration. For distributed tomography, each 
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node works on its own independent piece of the model. The overlapped pieces between 
adjacent nodes may produce different velocity estimates due to different input data on 
each node. To ensure the global smoothness of the velocity model, adjacent partitions 
need to exchange information from the overlapped areas. Global smoothing across 
5 computer nodes is illustrated in Figure 3. Data is exchanged between nodes 303 for 

smoothing every several iterations. As illustrated in Figure 3 control parameters are read 
from the deck file 301, and this information is shared across the nodes. Other parameters 
are read and computed 305. At 307 residual slowness is solved for and then smoothing is 
applied. The process at 307 iterates between solving for residual slowness and applying 
10 smoothing. Slowness values are output 309 after parameters have converged sufficiently. 



[0037] Global smoothing is necessary during velocity updating because we want the 
model difference between any two nodes to be small. Global smoothing enforces model 
consistency between nodes and global convergence of the solution. 

15 

[0038] Obviously, frequent exchanging of data between nodes increases network 
overhead. This overhead can be reduced by applying global smoothing every several 
conjugate steps. In addition, the data transmission between nodes can be reduced further 
by choosing a proper partition. 

20 

[0039] The general flow of this proposed distributed tomography is given in Figure 3 
with all the major features highlighted. Data communication is based on MPI. A Master 
node is responsible for distributing model and data to all other nodes. After each node 
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finishes its computation, it sends the updated local model to the master node. The master 
node assembles all the pieces and writes the result on to disk. 

[0040] Tomography on distributed PC clusters is an excellent way to solve large velocity 
5 updating projects with a reasonable computation time and relatively small cost. By 

partitioning a large model in to small pieces, a high density of input depth residuals 
becomes possible and thus highly accurate velocity inversions are possible. Issues 
associated with distributed computation need special treatment. There are areas 
overlapping between adjacent nodes which leads to model consistency among all the 
10 nodes. Area overlapping ensures the combination of inversion on each individual nodes 
is equivalent to the inversion of the entire problem on one node. Global smoothing 
provides an additional constraint for global convergence. Performance testing 
demonstrates that this strategy can reduce computation time significantly and increases 
throughput. 

15 

[0041] Figure 4 A illustrates a flowchart of a preferred embodiment of the method and 
system provided by the present invention. Prestack seismic imaging is performed 401 
using an initial velocity field and interpreted horizons. A residual depth difference is 
estimated 403 referenced to fixed offset and all horizons. Residual depth difference maps 
20 (p-maps) are computed 405 for each offset and each horizon. The residual depth 

difference maps are back projected 407 to determine slowness perturbation. The initial 
velocity model may be converted 409 to slowness and the estimated slowness is 
composited therewith to produce a new slowness volume. The new slowness volume is 
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converted 411 to a new input velocity volume for performing prestack seismic imaging. 
This process is repeated 413 until the slowness perturbation is negligible or reaches a 
predetermined threshold. Figure 4B illustrates an alternate embodiment of the invention 
where the threshold for slowness perturbation is determined prior to converting slowness 
to a new input velocity volume. 

[0042] The method and system of the present invention disclosed herein may be 
conveniently carried out by writing a computer program to carry out the steps described 
herein on a work station or other conventional digital computer system of a type normally 
used in the industry. The generation of such a program may be performed by those of 
ordinary skill in the art based on the processes described herein. Figure 5 illustrates a 
computer system comprising a central processing unit 1011, a display 1001, an input 
device 1021, and a plotter 1031. The computer program for carrying out the invention 
will normally reside on a storage media (not shown) associated with the central 
processing unit. The computer program may be transported on a CD-ROM or other 
storage media shown symbolically as storage medium 1041. 

[0043] The method and system of the present invention provides results that may be 
displayed or plotted with commercially available visualization software and computer 
peripherals. Such software and computer peripherals are well known to those of ordinary 
skill in the art. It should be appreciated that the results of the methods of the invention 
can be displayed, f lotted and/or stored in various formats. 
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[0044] While the invention has been described and illustrated herein by reference to 
embodiments in relation to the drawings attached hereto, various changes and further 
modifications, apart from those shown or suggested herein, may be made herein by those 
skilled in the art, without departing from the spirit of the invention, the scope of which is 
5 expressed in the following claims. 



COR-1060 



17 



